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Abstract 


An upwind MUSCL type implicit scheme for the three-dimensional Navier- Stokes equations is 
presented. Comparison between different approximate Ricmann solvers (Roe and Oslier) are 
performed and the influence of the reconstructions schemes on the accuracy of the solution as 
well as on the convergence of the method is studied. A new limiter is introduced in order to 
remove the problems usually associated with non-linear upwind schemes. The implementation 
of a “diagonal^ upwind implicit operator for the three-dimensional Navier-Stokes equations is 
also discussed. Finally the turbulence modeling is assessed. Good prediction of separated flows 
are demonstrated if a non-equilibrium turbulence model is used. 



1 Introduction 


Much efforts have been deployed during the last years to build efficient numerical methods for the 
solution of the equations of compressible viscous flows. Following the pioneer work of Godunov 
[16], efficient upwind methods were developed [42,36] which enable an almost perfect capture 
of stationnary shocks. With the introduction of essentially non-linear schemes, TVD [52,17,47] 
and ENO [19,18] schemes, the good shock capturing properties of the first-order upwind schemes 
were extended to higher order schemes. Acceleration techniques were also improved, reducing 
the cost of the computation of steady states. The diagonal dominance of the jacobian matrices 
obtained with upwind schemes was for instance exploited to construct relaxation based implicit 
methods [8,27,49], which may eventually be combined with multigrid methods [20,3,31]. Non- 
elliptic multigrid techniques were also introduced for the centered approximation of the Euler 
equations [21]. In this case, the damping of the high frequencies arises from a careful design 
Runge-Kutta time stepping scheme combined with an adapted non-linear artificial dissipation 
term [23]. In [28] the procedure was extended to the Navier-Stokes equations and the success of 
the method is clearly demonstrated in [40,54]. While the numerical techniques were drastically 
improved, in the meantime no much progress were reported in the modeling of the Reynolds 
stresses. Thus, even if the computation of the flowfield around a complete aircraft has become 
possible, [22,50], realistic simulation of separated flows around complex geometries remains 
unpractical. In this report a newly developed code, enabling the simulation of three-dimensional 
flows with reverse flow regions, on simple geometries (wings) is presented. The report is outline 
as follows; after a short description of the equations to be solved (§2), the numerical procedure 
is detailed: the spatial discretization in §3, the time integration and acceleration technique in 
§4 and the boundary conditions in §5. The turbulence models used are then presented (§6) and 
finally in §7 the results are discussed. 


2 Governing Equations 


The basic equations used as the physical model are the integral form of the mass-averaged 
Navier-Stokes equations — Reynolds equations. Since the dominant viscous effects for high- 
Reynolds number flows arise from the diffusion normal to the body surfaces, the thin-layer 
approximation is employed. The turbulent transport of momentum and energy due to the 
fluctuations of velocity and pressure is modeled with the eddy viscosity concept making the 
form of the Reynolds equations identical to the form of the Navier-Stokes equations. The 
equations to be solved are then 


i///v yrf ” + /iv Fii<i5 = /iv FvMs 


( 1 ) 


where 

U = ( p , pu, pv, pw, e) T 

is the vector of the conserved variables. The quantity V denotes an arbitrary control volume, 
dV and n correspond then respectively to the boundary surface and the outer normal of this 
volume. 
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In the cartesian frame the tensor F of the convection terms is 
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( 2 ) 


and the tensor Fy of the diffusion terms in the thin-layer approximation in the direction if, is 


Fy = 


'xy 


'xy 


' yy 


'yz 
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where e, represents tlie internal energy. The bulk viscosity A is evaluated using the Stokes 
hypothesis 

3 A 2 ft — 0 

and the molecular viscosity is determined from the Sutherlands law 


A Moo 


(£) 


f Too + 110.4 
T+ 110.4 ' 


(4) 


As a result of the eddy viscosity assumption, the equations (1) with the expression (3) for 
the tensor Fy correspond also to the Reynolds equations for turbulent flows if the molecular 
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viscosity is replaced by an “effective” viscosity // + /z ( and the quantity ^ by ^ The 

symbols fi t and Prt used in the previous expressions denote respectively the eddy viscosity and 
the turbulent Prandtl number. 


3 Numerical Method 

The equations (1) are solved with an implicit upwind method of the form 

(/ + At L)6U = -At R n 

f/ n+1 = U n + SU (5) 

where L is a spatial operator and R is the residual of the steady equations. No accuracy 
restriction will be imposed on the implicit operator as we are interested only in steady state 
solutions. 

3.1 Computation of the residual 

The computation of the residual is performed following the idea of Godunov [16] and gener- 
alized by Van Leer [53]. The procedure comprises of three stages. The first stage consists 
of the Reconstruction of the flowfield from its cell average values by piecewise polynomial 
approximations. In the second stage, the time Evolution (in the small, i.e. 0 < t < s) of the 
reconstructed flowfield is computed by solving, at each interface, a Riemann problem. The pro- 
cedure is therefore “upwind” as the wave propagation is taken into account when the Riemann 
problem is solved. In the final stage, the solution obtained by the resolution of the Riemann 
problems is Projected on the cells and new average values are computed. This step makes the 
overall procedure conservative. 


3.2 Reconstruction schemes 


Two reconstruction schemes were employed, a second-order ENO scheme of Harten and Osher 
[19] and a family of upwind biased schemes known as the k scheme. With the second-order 
ENO scheme, the smoothest parabola on the interval [x t -, x,+i] is first computed, 


P i+ i(*) = ?< + (x - x t ) 


<7.+i - <lx 


Ut+1 


+ (.t - ,r t+ i )t) 2 q 


with 


Dq t 


minmod(Dg,_i, Dqi) 

1 f Qi± 2 ~ g«+ 1 _ </«+! ~ <?. 1 

#t+ 2 “ J 


where the minmod operator is defined as 


minmod(.r. y) = 1 


sign(x) min(|x|, \y\ ) 

0 


if sign(x) = sign(y) 
otherwise. 
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The derivatives of P. , j.(x) and P._i(x) are then computed at a-, and the value of the slope of 
* ' 2 2 

the reconstructed field on x- i, x-. i is obtained from 

* 2 2 

. . rrfP, + i(a:.) dP t _i(z,)' 

S 1 q = minmod -2 , -j . 

(LX CLX 

The cell interface values are therefore, 

= V + ^ l 9 

tei-l = v-ji'i 


hi — x i + i x i-i- 

With the second-order ENO scheme, the computation of u” +1 involves the value of u n at 
the points i - 3 through i + 3 . This seven points stencil is required in order to distinguish 
between the variation of the flowfield near a shock from the variation of the flowfield near a 
smooth extremum. Five point schemes ( like the k scheme presented below ) cannot ensure, 
even if the limiter is designed such as it does not “clip” at an extremum, second-order accuracy 
in the vicinity of all smooth extrema. 

When the k scheme is used, the left and right values at the t + 1 cell interfaces are computed 
with an upwind biased interpolation, 


fci+i 


(1 + + (! ~ 

(1 + K)%t-1 +(1 -K)S<li+l 


Sq { .i = 


ft+i ~ Qi 
hi+i + hi 


In order to ensure an interpolation which does not increase the Total Variation of the initial 
distribution, the gradients 6q have to be limited. Following [6], the gradients were limited by 
replacing Sq by 6 with 

6+ = minmod [ ? >+1 , fr f* ~ r‘~ - 1 


6; = minmod 


h,+ 1 + h{ ’ h{ -|- h {- 1 

<?i - qi-i ^ gi+i ~ g» 
hi + hi-i ’ hi+i + hi 


where the value of b corresponds to the largest value for which the interpolated q- + i lies between 
q, and q,+i, 

. 3 — K 

0 — . 

1 — K 

The cell interfaces values are thus computed with 


= < 7 . + [(1 + k ) 6 ^ + (1 - k ) 6 ~ 

— 9i+l “ ^ 2 *- (1 + k)^i + (1 - 


*Li+\ 
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( 7 ) 


For k € [—1,1] the k scheme is second-order accurate and stable. The particular value 
k = 5 leads even to a tliird-order accurate scheme (on a uniform one-dimensional mesh). 
With the Chakravarthy-Osher limiter, the method is robust and when combined to the “third- 
order” scheme, fairly accurate results on coarse grids are obtained — see the result section. 
Unfortunatly, the Chakravarthy-Osher limiter has also the drawback of preventing the scheme 
to converge to the computer accuracy, hence a different limiter is proposed. The design of the 
new limiter was based on the following observations. 

• In order to have a good convergence rate, smoothness is the name of the game, thus all 
non-linear changes should be made as smooth as possible. 

• The only requirement in the design of a “monotone” shock capturing scheme, is to exclude 
the use of the gradients accross the shocks. 

The procedure developed consists therefore to locate the shocks and to modify, in a smooth 
way, the gradients in their vicinities. As the difference of the unlimited left and right pressure 
values at the cell faces 

x = \ P i+Ll - P i+LR\ 

is of the order of the truncation error in continuous flow regions and of the order of the shock 
strength in the vicinity of the shock, it can serve as a shock detector. The slopes are then 
limited according to 

s+ = h + ; ( 8 ) 

where $ is the restriction to the interval [0, 1] of the function a(x - e) 3 + 1. with 

a = -10 s 
£ = 1(T 7 . 

The coefficient £ represents the value over which the limiting is effective. As $'(£) = $"(£) = 0, 
the limiter is smooth and the scheme is not too much sensitive to the value of e. The coefficient 
a controls the strength of the limiting and was fixed by numerical experiments. 

An other way used to limit the influence of the gradient accross the shock is to combine 
a fully upwind reconstruction, k = — 1, with an upwind computation of the fluxes. As the 
gradients have to be limited only near the shocks, a natural procedure is then to blend, using 
the shock detector function $, the actual value of k with the fully upwind value k = — 1. 

k = + (1 - $)(-!). (9) 


In general, the downwind and upwind gradients should not be limited equally. For instance 

in the ideal shock case presented in the figure (1), only the gradient 6q - , i must be limited. In 

* ' 2 

order to distinctly limit the upwind and downwind gradients the ratios 


and 
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are introduced. If 7' < 1 (S + >• b ), respectively r + <C 1 (S 6 + ), only ^ + , respectively 

S~, must be limited. The restriction to the interval [0, 1] of the function 


*> = ; <’'«> = °» ( 10 > 

is therefore proposed to compute the limited gradients 6 ± , 

6* = [$ + (\-*)e(r ± )}6 ± . ( 11 ) 

The value 0(1) has to be set to zero in order to limit equally and S~ when the downwind 
and upwind gradients are identicals. This choice fixes the value of b. The coefficient c was 
determined by imposing 0(0.5) = 0.5. In all the computations done with this new limiter, the 
value of e was set to 0.003. 

In summary, with the “smooth limiter” just described, the values q i+ i at the cell faces are 

computed using the relation (7) in which the gradients were obtained through (11) and 
where the value k was replaced by the blending k given in (9). 


3.3 Resolution of the Riemann problem 

In the original work of Godunov, the Riemann problems were solved exactly, but it has been 
proven that similar flowfields could be obtained at a lower cost if the Riemann problems were 
solved only approximately. Two of such approximate Riemann solvers are those proposed 
by Roe [42] and Osher [36]. Both rely on the wave decomposition. In the Roe scheme the 
decomposition is made by assuming a locally constant state U for which 


Ur-Ul = Hk^k^k 
F(Ur)-F(U l ) = E k^khotk 

where 

Afc is an eigenvalue of the jacobian matrix ff(ff) 

f*. is a right eigenvector of §p(U) 

f ~*k — h(UR — Ui) with Ik a. left eigenvector of | jj{U) 


In [44], Roe and Pike proved that the state U for which the equations (12) are satisfied exist 
and is unique. It can be computed using the special averaging introduced by Roe in [42], 


P 

u 

v 

w 

H 


<xpn 

1+a 

ovl±VJ2 

1+a 

atcr+TUp 

1+a 

1+a 


( 13 ) 
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where 


“ = vf 

H = ~l + ^(w 2 + v 2 + to 2 ). 

The Osher scheme is in some sense, a non linear extension of the Roe scheme. In the case 
of a non-linear hyperbolic system, the wave decomposition can be written as, 

Ur — Ur = Ylk Ir k r k . . 

F(Ur)-F(U l ) = Ekfr k hnds k 

where the T k curves are defined by 

dU 

drr n - 

Using the constancy of the Riemann invariants along the T k curves, explicit expressions for the 
integrals can be obtained [36] [29], 

Jr^rids! = F(U t+ x) - F(Ui) 

/r 2 A 2 r 2 ds 2 = F(U i+ z)-F(U i+k ) (15) 

fr^radss = F(U i+1 )-F(U i+i ) 

with 
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ft»,-+n,-+i 
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2 a,. t-i-afl; 
7— 1 1+a 


1 : 


w i+i 




r— 1 | q t4-i^ a > 

2 1+a * l+o 

'±Lil* 

7 Pi 

' a 2 

,_!±i 

*+3 7 


an,+n, + i , 2 rt,+i -oa,: 

1+a ’ 7— 1 1+a 


P variant 




a i+l 




P l+ i 


7-1 Mi-tK-n I a, ■+!+«,' 

2 1+a ' 1+a 


W+iV - 1 


O variant 


a !i±i 

fttiT 


w «+l 

= 

v i+l 

= v i+ 1 

w i + 1 

= Wi+1 

a i+'i 

= aa i+ 

Pi+I 

- P J±i 

~ a 2 

p iH 

= p i+i 


P and 0 variant. 


Intermediate states between U { , U i+ 1 and fA + 2 , J 7 l+ i have also to be introduced in order to 

account for a possible change of the sign of A = u + a or A = u - a. These intermediate states 
are 
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> P variant 



n *+* 


= 2+f(«i + 73r) 
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W U 

°ui 


PL i 


P* 


= Wi 


u ‘+i 

£!*±i 
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p Ui 
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/ 0 variant 
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> P variant 
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O variant. 


( 19 ) 


Both with the Osher and Roe scheme the flux at the interface is computed using, 


with the AF 4 e fincd as 


and where 


F i+i = f<Pl) + 

k 



or 


A-|A| 
A --y-. 


( 20 ) 


( 21 ) 

( 22 ) 


Because in the equation (21) the integration is made only on the negative part of A, the 
integrand is not an exact differential and the values of the fluxes depend on the ordering in the 
integration path (JIV For the Osher P variant, the integration path is r u -o , T u , T u+a and for 
the 0 variant, it is r„+ a , T u , r u _ a . A closer inspection of the integrals shows that if unlikely 
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cases are discarded - supersonic flows in opposite directions on each side of the interface for 
instance — the computation of the flux at the interface requires at most three evaluations of 
fluxes for the P variant while it requires four evaluations for the 0 variant. 

The above expressions are useful for computations in a cartesian mesh ( the splitting for the 
Osher scheme was made in the i direction ). For general meshes, the splitting is made in the 
direction normal to the cell face. Such a choice is computationally convenient but is also quite 
arbitrary; the splitting direction depends on the mesh rather than on the flow properties. In 
practice however, the procedure appeared to perform well and even if different grid independant 
splitting formulations have been known for several years [13,43,14,38], their efficiency are yet not 
fully established. With the “grid splitting” procedure, the extension to generalized coordinates 
of the approximate Riemann solvers is simple, the cartesian velocities u , v and w have just to 
be replaced by the contrava riant velocities u , v and w. 

u = n x u + n y v + n z w 

V = l x u + l y v + l z w (23) 

w = m x u + m y v + m z w 

where n is the vector normal to the cell interface and l , m are the tangential vectors. In 
the two-dimensional case, the tangential vector is uniquely defined by the knowledge of the 
normal vector, in the^three-dimensional case however, the normal vector defines only a family 
of tangential vectors / and m. A precise definition of these tangential vectors can yet be avoided, 
by noticing that only the linear combination of the tangential velocities 

fi* = l x v + rrtxW 

Sly — l y v -f rriyW 
Sl z = l z v + m z w 

are needed, and that these linear combinations can also be written as 

= u — n x u 

Sly = v-Uyii (24) 

Sl z = w — n z ti. 

Thus in practice, the tangential velocities v , w arc replaced in the computations by the velocities 
Slj? , Sly , and Sl z • 


3.4 Computation of the viscous terms 


In a finite volume approach, the computation of the viscous terms requires the evaluation of 
first-order derivatives at the cell faces. The computation of such gradients may be obtained 
using the Gauss theorem [37,35] but in this work a simpler procedure has been employed and 
the gradients were computed with 


dq. 


(?«+ 1 - < 7 <)£,+ 1 


"H-i 


+ V, 


(25) 


where is the surface of the interface and V t , V;+i are the volumes of the cells on both 
sides of the interface. 
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4 Implicit Operator 


The implicit operator used to accelerate the convergence to the steady state is an extension 
to the three-dimensional case of the operator presented by Coakley [10]. It is derived from a 
backward Euler implicit integration of the equation (1). Instead of solving this operator by 
relaxation techniques as done in two dimensions by [27,49,30] and partially in three dimen- 
sions by [7], the three-dimensional implicit operator is dimensionally splitted [5] leading to the 
resolution of three one-dimensional operators of the form 

(/ + = R HS (26) 

where A is the jacobian matrix of the Euler flux in the direction £ 

- OF OF 0C OH 

A ~ d£ ~ ou n x + ou Hy + ou Hz 

with 

F = Ff; G = F j ; II = Fk 
and M is the matrix of the viscous terms. 

If the spatial derivatives of (26) are discretized with a three points stencil, the resolution of 
the system (26) necessitates the inversion of a bloc 5x5 tridiagonal matrix. Following Cha.ussee 
and Pulliam [9] the system can be “diagonalized” leading to the inversion of scalar tridiagonal 
systems. This is accomplished by replacing the matrix M by its spectral radius <7 J, using the 
decomposition, A = T -1 AT and by taking the product T _1 A out of the spatial derivative. 
Introducing the characteristic variables SU = TSU and multiplying the equation (26) by the 
matrix T, the system to be solved is now, 

(I + At -At a / J j;)6U = RHS. (27) 

The Coakley scheme is finally obtained if a first-order upwind approximation is used for 
the spatial derivatives of the convective terms and a centered approximation for the diffusion 
terms. It should be noted that whereas the matrix A involves only the components of the normal 
vector n, the matrices T~ l , T and thereafter SU involve also the tangential vectors / , m. In 
consequence, while the original system requires only the knowledge of the normal vector n, the 
“diagonalized” system necessitates also the definition of the tangential vectors Z , m. As the 
normal vector ft was calculated by taking the cross product t\ x f 2 fig. (2), it seems natural to 
choose 


f= -4r- 
\M 

(28) 

m = n x /. 

(29) 


Such a procedure however produces spurious crossflows for purely two-dimensional flows. In 
the simple case of a two-dimensional turbulent flow on a flate plate, the vectors J, k can be 
chosen as the tangential vectors if the mesh is cartesian, fig. (3). The contravariant velocities 
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are then identical to the cartesian velocities, 


Am = Am 
Am = Aw 
Atw = Au> (= 0) 

and the implicit system for the tangential velocities can be written as 

C[6t>] = [Aw] 

£[£tw] = [Atw] 


( 30 ) 


rp 

[H = [$Wj i , ivj—i , 6vj , ^Wj.|.i , • • • , Svj ] . 

If the matrix £ is not singular, no crossflow can be generated by the implicit operator. 

With the procedure proposed for the computation of the tangential vectors l and m, equa- 
tions (28,29), the vectors will differ from /, k, there form will be 


h = ctjf + Pjk 

m j = IjJ "h ^3 ^ • 

If these vectors are used for the definition of the contravariant velocities, the “diagonal” implicit 
system can be written as 

C[6v\ = [a] [Aw] + [/3][Aw] 

C[f>w\ - [7] [Aw] + [<$][Aw] 

thus 

IM = r-^alfAw] + £ _1 [/?][Aw] ^ 

[tftw] = £ -1 [7 ][Aw] + £ _ 1 [£][A«;]. 

The solution of the implicit operator will be independant to the choice of the tangential vectors 
if the tangential velocities satisfy, 


[H = [apw] + [/?][<H ^ 

[£tw] = [t][^w] -f [£][£«;]. 

Unfortunatly the relations (32) are not always fulfilled as it can be seen by replacing in (31) 
[Aw] by £[£w] and [Aiw] by £[Atw], 

m = c-'wcM + c-'imn (33) 

[Sib] = £- 1 [7]-C(H + £-* [$]£[$«>]• 
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The expression (33) can be identical to (32) if and only if the matrix C commute with the 
diagonal matrices [a] , [/?] , [7] , [$]. For general matrices £, this necessitates that 


[а] -- a I 

m = (31 

[7] = 7 / 

[б] = SI. 


(34) 


In the case of the two-dimensional turbulent flow on the fiate plate, because of the stretching 
of the mesh in the boundary layer, the relations (34) are far from being satisfied and as a 
consequence, a strong crossflow is generated, 6w ^ 0. The “natural” procedure has therefore 
to be rejected. 

The proposed method for the computation of the tangential velocities, follows the procedure 
used at the explicit step. The computation of the increments of the tangential velocity vectors 
6v and 6w are replaced by the computation of linear combinations of them i.e. 


6$^ = ( l x bi) + m. x 6w) — 6u - n x 6u 

fifty — (lyfv + rriySib) = 6v — n y 6u (35) 

fft z = (l z 6v + m z 6w) = 6w — n z 6a. 


Performing some linear combinations on the equations for 6v and Sib, the equations (36) can 
be formed, 

6% + uA t = AQo ; 6 = x, y, z. (36) 

Finally the equations for 6 ft# are obtained if the coordinate components l#. m # are taken inside 
the spatial derivatives — this approximation is of the same order as the approximation used in 
the “diagonalization” — leading to 

Q 

(I + uAi — )6fte ~ Aft# ; 9 = x, y, z. (37) 

Thus instead of solving the two equations for 6v and fiw, the three equations for 6ft x , and 
Sft z will be solved. Of course, the three velocity components Sft# cannot be independant, they 
should satisfy 

n x bft x + Tiybfty + n z 6ft z ~ 0. (38) 

Because of the simplification made for the transformation of (36), the equation (38) is generally 
not verified ; the tangential vector 6ft does not lie in the tangential plane. Therefore it is 
replaced by its projection on the tangential plane 6ft\ 


6ft f Q = 6ft $ — n$6 


(39) 


with 


6 = n x Siil x + HySCly + n z 6 Q, z . 


By construction no crossflow can be generated when this procedure is used for the computation 
of the implicit tangential incremental velocities. Numerical experiments have also shown that 
for purely two-dimensional flows, the three-dimensional operator leads to almost the same 
convergence rate than the corresponding two-dimensional operator. 
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5 Boundary Conditions 


The use of upwind schemes simplifies significantly the computation at the boundaries as solely 
the “physical” boundary conditions are needed. For instance, if only simple waves are assumed 
at the boundaries, the Osher P variant scheme can be used at the boundaries by replacing 
the information contained in the waves outside of the computational domain by the natural 
physical boundary conditions [15] : static pressure at an outflow boundary; rest temperature, 
rest pressure and flow angle at an inflow; no slip velocity, zero normal pressure gradient and 

adiabaticity at a solid wall For lifting airfoils, the freestream conditions are modified by 

taking into account the circulation around the airfoil as discussed in [48]. The undisturbed flow 
conditions are used for flows around finite wings as the disturbances created by the wing decay 
more quickly in the three-dimensional case than in the two-dimensional one. If sufficiently far 
from the wing, they behave like those generated by a point singularity rather than by a line 
singularity. 

At the implicit step all the boundary conditions are treated implicitly. For the characteristic 
variables SQ X , Sil y , SQ Z and a 2 Sp — SP, if the index of the boundary is denoted by N + 1, the 
equations at all type of boundaries can be cast into the form 


= diSUpt + d 2 SUw-i. 


(40) 


with 


di = d 2 = 0 for a solid wall, a farfield boundary 

di = | , d 2 = — ^ for an extrapolated boundary condition 


d\ = 1 or 0 , d 2 — 0 for a symmetry condition 


If treated implicitly, the symmetry and wall boundary conditions introduce some coupling 
between the variables SR + = SP + paSu and SR~ = SP - paSu. The implicit equations can 
therefore not completely be decoupled and a 2x2 block tridiagonal system has to be solved for 
these characteristic variables. The boundary conditions associated with this 2x2 system are the 
equivalent of (40) 


6R+ 


dt 


SR- 

N+l 

d t 

d 2 j 


with 

df, 2 ,3 = 0 

df = 4 = i, 4 = 0 

df ~ d 2 — f ' d 2 — d^ = 0 , f/3 = — ^ 


SR+ 

+ 

di 

0 


6R+ 

SR- 

N 

0 

1 


SR- 


for a farfield 

at a wall and a symmetry plane 

for an extrapolated boundary condition 


As the same functional form (40) or (41) is used at all type of boundaries, the computation of 
implicit boundary conditions does not impaired the vectorizability of the implicit step. 
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6 Turbulence Modeling 


Two simple mixing length turbulence models were examined in this work. The first model 
considered is the standard Baldwin- Lomax model [4] 

vu = pK 2 rj 2 D 2 \u\] K = 0.4 

v to = 1.6F to Fkleb(»;) (42) 

v t = mm(v ti ,v to ) 


where 
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= 

0.0168 [l -f 5.5(^^-) 6 ] _1 

L *?max ' J 


The definition of u T was chosen slightly different from the wall shear stress as used in the original 
Baldwin-Lomax turbulence model. This modification was introduced in order to prevent from 
the computation of a vanishing eddy viscosity in a section emanating from a saddle separation 
point (|w|w = 0). 

It is well known that equilibrium models such as the Baldwin-Lomax model, are not suited 
for separated flows for which the diffusion and the convection of the turbulence are no more 
negligible and introduce some imbalance between the production and the dissipation rate of 
turbulence. While retaining the eddy viscosity assumption these non equilibrium effects can 
be taken into account by two-equation models, K — e, K — ui, but it seems that despite their 
“universality” the two-equation models does not improve significantly the agreement between 
the computed results and the experimental data for separated flows [11]. A less ambitious 
approach is to modify two-layers mixing length models in order to extend their successes to 
separated flows. Such an approach was taken by Johnson and King [25] and the model they 
derived, appeared to be adequate for the computation of separated flows on airfoils and wings 
[11] [2]. 

The idea behind the Johnson-King model is (i) to scale the turbulent velocity to the square 
root of the maximum Reynolds shear stress rather than to a length scale wall vorticity product; 
(ii) to compute the maximum Reynolds shear stress by solving a differential equation in which 
non-equilibrium effects are taken into account. As the level of the turbulent shear stress is then 
determined by the differential equation, the Johnson-King model in contrary to standard mixing 
length models, neither depends only on local mean flow gradients nor assumed a turbulence in 
equilibrium. The eddy viscosity distribution in the inner layer used with the Johnson-King 
model is then 

i'ti = Kr)D 2 UM (43) 
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where 


*M = (“)* 

D = 1 - exp(-^i) 

+ _ nmax(uM,«r) 

Here and below the index a/ indicates the location where the Reynolds shear stress is maximum. 

In the original formulation of Johnson and King the outer eddy viscosity layer was based on 
the Cebeci-Smith distribution. This formulation was well suited in the boundary-layer context 
used by Johnson and King for the derivation of their model, but with Navier-Stokes codes the 
Baldwin-Lomax formulation is more convenient. The outer eddy viscosity layer is therefore 
calculated using 

i/ io = irl.bFtuFk (44) 


with 


6 — 1.9 T) max 

Fk = 0.0168 [l + 5.5(|) 6 ] _1 . 


The coefficient a is introduced to force the value of the maximum shear stress tm — to 

match the value t\j obtained by the resolution of the differential equation (48). This coefficient 
can be computed by solving the equation 


tm{o) -tm = 0 

with a Newton method, or with a procedure proposed by Abid [1], 


a t+At = g-t TM- 




= 1. 


(45) 


(46) 


Knowing the values of Vti and t/ io , the actual value of the turbulent eddy viscosity u t is 
computed with 

vt = /'(<,( 1 - exp(— — )). (47) 

Vto 

As stated, the level of the turbulent shear stress is obtained through the resolution of a differ- 
ential equation, see [25] for its derivation. In the three-dimensional case, this equation is 


. d x , . d x 0.125 
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\/ (1 — cri )+ 


(48) 


with 

L\j = min(KriM, 0.225 K6). 


The left hand side of (48) represents the convection of the turbulent shear stress. The diffusion 
of the turbulent shear stress being modeled by the last term of the right hand side, the remaining 
term corresponds then to the imbalance between the production and dissipation of turbulence. 
This term is consequently approximated by the difference between the actual shear stress and 
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the shear stress that would have been obtained if the turbulence was in equilibrium. This shear 
stress, r eq , is computed at t) — i]m using 


^eq — /J<eq|^| 

i'teq = "toeq (l - eX P (-^)) 

Vtit q = Kr)D 2 u eq 

Utoe q = l.GFtyFk 


(49) 


with 


D = 1 - exp(— 


>?«eq 

llu 


) 


We q = maX^teqM^M). 


The convection terms of (48) are approximate by first-order upwind differences. The equation 
is then solved with a Point Alternate Symmetric Gauss-Seidel Relaxation with one relaxation 
performed in each direction. An alternative approach used in [39] and [2], is to add a time 
dependent term and to solve the equation in the same way as the Navier-Stokes equations (1). 


7 Results and Discussions 

The newly developed code was first validated on two-dimensional airfoils and the influence of 
the different choices of Riemann solvers and reconstruction schemes was studied. For instance 
the influence of the Riemann solver on the solution was studied on the CAST10 airfoil at 
Mqo = -767, a = 2.159, Re = 1 X 10 7 . Under these conditions the shock is strong enough to 
induce a separation of the boundary layer. As it can be seen from figures (4-6), the results 
obtained with the Roe and Osher schemes are almost identical for this sensitive test case, even 
on coarse meshes. In particular the rapid acceleration at the leading edge is well predicted. 
The difference in the shock location between the experiments of [34] and the computations, can 
be due to the extreme sensitivity of this airfoil to wind tunnel side-wall effects [33], or are due 
to the Baldwin-Lomax turbulence model employed in the computations, model which is not 
adequate for separated flows. A clear improvement in the shock location was obtained with 
the Johnson-King model, fig.(7). Not only the computed pressure distributions on the airfoil, 
obtained with the different approximate Riemann solvers are the same, but the whole flowfields 
look alike fig.(8-10). 

The numerical procedure shows a larger sensitivity to the reconstruction schemes than to 
the Riemann solvers. On both the CAST10 and the RAE2822 airfoils, the inherent numerical 
viscosity of the second-order ENO scheme induces some significant differences, at the leading 
edge and on the upper surface of the airfoils, between the computed pressure distribution and 
the experiments of Mineck [34] and Cook et al. [12], fig. (11, 12). A better agreement with the 
experimental data and a solution closer to the fine grid solution was obtained if the third-order 
k scheme with the Chakravarthy- Osher limiter is used, fig.(4,13). It can be noticed however, 
that the pressure distribution on the lower surface of the airfoil is better predicted with the 
uniformly second-order scheme than with the TVD scheme which becomes a first-order scheme 
at extrema. Comparing the result obtained with the unlimited third-order scheme, fig.(14), and 
with the TVD scheme, it is clear that the limiter is active not only near the shock but also on 
the lower surface and near the leading edge. 
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In order to establish the optimal method, the computing time and the convergence histories 
of the different approximate Riemann solvers have to be compared. From table (1) it is clear 
that the Roe scheme is less expensive per iteration than the Osher scheme. The convergence 
histories have also to be examined as it has been found [41], that the Roe scheme induces 
some oscillations for slowly moving shocks whereas the Osher scheme does not generate such 
oscillations. As the shock approaches its final location, a slow moving shock is likely to appear. 
Therefore, if the oscillations are not damped efficiently, the convergence may slow down. This 
behavior has not been observed and even the convergence with the Roe scheme was found to 
be slightly better, fig. (15-17). Thus, as the third-order scheme with the Chakravarthy-Osher 
limiter, not only is more accurate than the second-order ENO scheme but is also less expensive, 
table(2), the optimal choice is clearly to combine the Roe scheme with the third-order k scheme. 

As already observed [32], the limiters have usually an adverse effect on the convergence rates 
of the scheme. This undesirable property can be illustrated by comparing on the RAE2822 airfoil 
at Moo — -676, a = 1.93, Re = 5.7 x 10 6 (case 1), the convergence history of the limited scheme 
with that of the unlimited one, fig.(18-19). The results indicate that even if the flow does not 
present any shock, the limiter is active and the non-smoothness of the limiting function creates 
some strong non-linearities which are responsible of the deterioration of the convergence history. 
It has been shown [55], that even with a differentiable limiter $(r), where r is the ratio of the 
successive gradients, the problem remains. The reason lies in the behavior of r which is like a 
random function in the farfield, where the gradients are very smalls. The cure in this case is 
to use a cutoff value under which the gradients are not limited, as described in [51]. It can be 
noted, that although only a two order of magnitude decay of the residual is obtained with the 
limiter turned on, the solution itself is almost identical to the unlimited one and agrees well 
with the experimental results, fig.(20-21). Therefore a convergence to the level of the truncation 
error can be suspected. It was found by numerical experiments, that a faster convergence is 
obtained when higher CFL numbers are combined with an under-relaxation of the increments. 
This behavior is illustrated by comparing the figures (19) and (22). In both cases near optimal 
CFL numbers were used. 

If the solution contains a shock, the unlimited third-order scheme does not lead to a smooth 
solution, fig.(23). The non-smoothness of the solution has also an adverse effect on the con- 
vergence, fig.(24). In such case, a smooth solution without limiter can be obtained if the fully 
upwind scheme (k = -1) is used, fig.(25). This value of k leads to a much better convergence 
history as well, fig.(26). With the limited scheme, a limit cycle was again obtained, figures 
(27-29). It is interesting to note that the level of the residual increases with increasing CFL 
numbers. With higher CFL numbers, the lift reaches its mean value more rapidly but noticeable 
fluctuations around the mean value become more and more apparent. When highly non-linear 
schemes are used to compute steady solutions, one should therefore ensure that the residual has 
effectively been reduced below the level of the truncation error. 

The convergence problems found with the Chakravarthy-Osher limiter can be elimitated 
on the RAE2822 airfoil, if the limiter is replaced by the “smooth limiter” described above. A 
convergence to the machine zero is then obtained on coarse and fine grids, fig.(30). It can be 
noted that the convergence rates of the scheme with the “smooth limiter” and of the unlimited 
scheme are alike; they even are almost independant to the grid density, fig.(30, 31). If the 
convergence properties of the “smooth limiter” and of the unlimited scheme are similar, the 
solutions calculated with the “smooth limiter” are more accurate than the solutions obtained 
with the unlimited scheme. Examining the figures (32, 33), it is clear that the solutions com- 
puted with the “smooth limiter” do not present any oscillation. Comparing on the RAE2822 
airfoil the pressure distributions on the coarse grid, it appears that the accuracy of the scheme 
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with the “smooth limiter” lies between the accuracy of the unlimited second-order scheme and 
the accuracy of the third-order scheme — with and without the Chakravarthy-Osher limiter — 
fig. (32). The pressure distribution with the “smooth limiter” for instance presents an overshoot 
at the leading edge whereas the third-order scheme does not show any overshoot. This over- 
shoot is however more pronounced with the fully upwind scheme. Similar results are obtained 
for the pressure distribution on the upper surface; the pressure level is better predicted with 
the “smooth limiter” than with the fully upwind scheme but the shock location and the pres- 
sure peak near the leading edge, are not as well predicted with the “smooth limiter” as with 
the third-order scheme. These results could have been expected as the value k is a blending 
of the value « = -1 (fully upwind) and k = | (third-order). It is interesting to note that in 
contrast to the Chakravarthy-Osher limiter, the “smooth limiter” does not smear the extremum 
on the under surface of the airfoil. Thus, even if the goal consisting of limiting the gradients 
only in the vicinity of the shocks, has not been fully achieved, significant improvements were 
obtained. The excessive dissipation introduced at the leading edge is a consequence of the lack 
of resolution in this region. As the mesh is refined, the results become more and more alike, 
fig.(-33). A closer inspection on the skin friction coefficient shows however, that the fully upwind 
scheme introduces more dissipation at the trailing edge than both the third-order scheme with 
the Chakravarthy-Osher limiter and the “smooth limiter”. On the skin friction coefficient, the 
shock appeared also to be better captured with the “smooth limiter” than with the fully unlim- 
ited upwind scheme. From these experiments one can consider that the third-order scheme with 
the Chakravarthy-Osher limiter is the most accurate scheme used, but the k scheme with the 
“smooth limiter” while maintaining most of the accuracy of the third-order scheme, improves 
drastically the convergence of the scheme. 

The “smooth limiter” was also applied to the computation of the flowfield around the same 
RAE2822 airfoil but with different flow conditions, = .75, a = 2.81, Re=6.5 x 10 6 (case 
10). Under these conditions a shock induced separation forms on the upper surface of the airfoil. 
With the Baldwin-Lomax turbulence model, a wrong shock location is again obtained fig.(34) 
but no oscillations are created and the “smooth limiter” does not produce too much dissipation 
at the trailing edge of the airfoil, fig.(34). In this case also, a machine accuracy convergence 
was reached, fig.(35). On the ONERA M6 wing M 0 0 = .84, a = 3.06, Re = 11 x 10 6 , the 
“smooth limiter” improves the convergence of the method when compared to the Chakravarthy- 
Osher limiter, but the convergence rate is not as good as that of the unlimited fully upwind 
scheme, fig.(36). Examining the solution, fig.(37), it is apparent that the results calculated 
with the “smooth limiter” are in a large part of the wing, close to the results obtained with the 
Chakravarthy-Osher limiter. While resolving slightly better the leading edge pressure peak, the 
unlimited scheme smears much more the shock. Near the tip of the wing, ^ > 0.9, the “smooth 
limiter” introduces an additional extremum. This new extremum is probably responsible of the 
slowdown of the convergence as compared to the unlimited scheme. It is the author belief, 
that if the limiter is designed correctly, the convergence with the limited scheme should be as 
good as the convergence obtained with the unlimited scheme. The behavior of the scheme with 
the “smooth limiter” on the coarse grid indicates that something was not done properly. This 
indication has been confirmed by the difficulties encounted on a finer mesh with the “smooth 
limiter”. In conclusion, the “smooth limiter” improves significantly the convergence of the 
method on the RAE2822 airfoil, but it is still not the “Ultimate” limiter and modifications 
have to be made in order to obtain a broader range of its applicability. 

The influence of the turbulence modeling has finally been studied. For attached flows on 
airfoils, we have seen fig.(18,33) that a good agreement with the experimental data can be 
obtained with the Baldwin-Lomax turbulence model. The same is true for the attached flow on 
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the ONERA M6 wing = .84, a = 3.06, Re = 11 x 10 6 , if the grid is fine enough, fig. (38). 
For separated flows, (RAE2822 case 10), the Baldwin-Lomax model predicts a shock location 
which is far downstream of the location observed experimentally, fig.(39). No improvements 
are found by refining the mesh. It can be noted once again, that the third-order scheme 
combined with the Chakravarthy-Osher limiter gives fairly accurate results on coarse grids. If 
the equilibrium model does not predict the proper shock location, with the non-equilibrium 
model of Johnson and King, the correct shock location is found, fig.(40). Comparing the grid 
refinement study made with the Baldwin-Lomax model and with the Johnson-King model, 
the greater sensitivity of the Johnson-King model clearly appears. A noticeable difference in 
the shock location is observed whether a 161x33 or a 257x65 mesh is used. Futhermore, the 
Johnson-King model is particulary sensitive to the initial level of the maximum shear stress. For 
instance the “solution” (unsteady) shown in figure (41) was obtained in the following manner : 
the Reynolds equations with the Baldwin-Lomax model were first solved on a 81x17 mesh; this 
solution was then interpolated on the 161x33 mesh and serves as the initial condition for this 
finer mesh; fifty iterations were then performed on the 161x33 mesh with the Baldwin-Lomax 
model before the original Johnson-King model was turned on. If more iterations had been 
performed on the 161x33 mesh with the Baldwin-Lomax model, steady but wrong solutions 
would have been obtained as well. The solution shown in figure (40a) was computed by taking 
as the initial condition, the converged solution obtained with the Baldwin-Lomax model on 
the 161x33 mesh. In order to enforce a unique solution — independant of the initial condition 
— the velocity scale used in the inner layer of the equilibrium eddy viscosity, has to be 
replaced. Instead of using 

n e q = I'teqM 

as in the original Johnson-King model, 

U eq = max(l/* eq M, l/\u\) 

must be employed. This fix was yet not sufficient in the computation of the ONERA M6 wing 
shown below. The level of the starting maximun shear stress was in this case still to low and 
in consequence, the shock location was moving upstream without any bound. The fix was then 
to replace in the computation of t'tieq, the inner layer formulation of Johnson and King by the 
Baldwin and Lomax formulation, 


I/f,eq = K 2 T] 2 D 2 \U)\. 

With the Baldwin-Lomax formulation for v t ie q , a better solution was found on the RAE2822 
(case 10) airfoil on a coarse grid, fig.( 42), but on a finer grid, the pressure recovery was 
unfortunately not predicted as well as with the “original” formulation; a pressure bump is found, 
fig.(43). This bump has also been observed by Radespiel [39] and Swanson [46], with a Jameson 
type scheme and with the “original” formulation of Johnson and King. The extreme sensitivity 
of the Johnson-King model shows up also on the convergence of the method. From figure (44), 
it is apparent that the lift converges with difficulties, the clear decay of the oscillations indicates 
nevertheless that if enough iterations are performed, a steady solution can be expected. It should 
also be pointed that the computation of a proposed by Abid, equation (46), introduces some 
time dependency in the solution. The steady state will slightly depend on the time integration 
path and depending whether the turbulence quantities were updated at every iterations or only 
at every five iterations, different convergence histories were found. 

As already noticed by Coakley [11], we also observed that whereas the Johnson-King model 
improves significantly the prediction of separated flows, attached flows are not as well resolved 
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with the non-equilibrium model as with the Baldwin-Lomax model, fig. (45). The reason again 
lies in the inner layer eddy viscosity formulation of Johnson and King. Recently, Johnson and 
Coakley [24] proposed a new formulation which consists of a non-linear blending of the Johnson- 
King and Baldwin-Lomax formulations. This new law seems to ameliorate the results but the 
added non-linearity increases also the convergence problems. 

The behaviors of the Johnson-King model observed in two dimensions were confirmed in 
three dimensions. The test case used is the ONERA M6 wing, where comparisons with both 
detailed experiments [45] and numerical results [2] [39] are possible. For an angle of attack 
of o = 6.06, a large separation region forms on the upper surface of the wing, fig. (46). In 
this case the Baldwin-Lomax model is again inadequate to capture the main feature of the 
flow. The pressure plateau after the interaction of the two shocks — the shock emanating from 
the leading edge and the normal shock, fig.(47) — is for instance not captured at all and this 
disagreement with the experimental data is not due to a lack of resolution, fig.(48). As pressure 
plateau regions behind a shock wave are usually the result of a large reverse flow region, a non- 
equilibrium model must be used in order to obtain a good resolution of the pressure distribution. 
With the non-equilibrium model of the Johnson and King, it is clear that a good representation 
of this pressure distribution is obtained, fig. (49) and the expected large reverse flow region 
is effectively found, fig. (50). These results are the consequence of the lower values of eddy 
viscosity predicted by the Johnson-King model in adverse pressure gradient regions. The lower 
values of eddy viscosity also induce an upstream movement of the shock location, fig. (51). No 
experimental visualizations of the wall streamlines were available on the ONERA M6 wing, but 
the computed, mushroom type, wall streamlines are in good agreement with the results obtained 
by Abid, Vatsa et al. [2]. The mushroom type structures were also observed experimentally on 
other wings [26]. If the Baldwin-Lomax formulation for ^ <ieq is used, computation on coarse 
grids with the Johnson-King model were possible and satisfactory results found, fig. (52-54); 
even the wall streamlines patern was qualitatively well predicted, fig. (54). For an attached 
flow case, a = 3.06, the Johnson-King model predicts a shock location which is as in the two- 
dimensional case, slightly upstream of the position obtained with the Baldwin-Lomax model. 
The difference however seems to be less pronounced in the three-dimensional case than in the 
two-dimensional one, fig.(55,56). 


8 Conclusion 

Two and three-dimensional computations have been presented and differences between several 
upwind schemes discussed. It has been found that the differences due to the upwind schemes are 
negligible if the three waves existing in the Euler equations are taken into account by a Riemann 
solver. The reconstruction scheme was found to have more influence on the accuracy of the 
solution. Accurate solutions on coarse meshes were obtained with a third-order upwind biased 
scheme. Such scheme requires the use of some limiter in order to compute smooth solutions 
with shocks. If the use of limiters leads to robust schemes, it has the drawback of preventing 
a convergence to the machine accuracy. It was proven that the use of a “smooth limiter” 
can drastically improve the convergence of the method. The “smooth limiter” proposed, while 
leading to good convergence rates on the RAE2822 airfoil, presents still some defects on the 
ONER A M6 wing. Thus modifications on the limiter need to be done to enlarge its applicability. 

Equilibrium and non-equilibrium mixing length type turbulence model were tested on two 
and three-dimensional configurations. While the equilibrium model was found to give accurate 
results for attached boundary layer type flows, it was also proven to be inadequate to pre- 
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diet correctly separated flows. Inversly, the solutions obtained with a modified version of the 
Johnson-King model agree well with experimental data for separated flows, but are less accurate 
than the solutions computed with the Baldwin-Lomax model for attached flows. Futhermore, 
an excessive sensitivity of the non-equilibrium model was experienced. Therefore some changes 
in the non-equilibrium model in order to improve the accuracy of the solutions for attached flows 
and to reduce the sensitivity of the model, have to be made before it can be used routinely in 
engineering computations. 
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Figure 3: Evolution of the tangential vector m in the boundary layer 



Figure 4: 

Pressure distribution 
Osher 0 scheme, k = ^ 
Chakravarthy- Osher limiter 
81x17 mesh (48 points on the airfoil) 



Pressure distribution 
Osher P scheme, k = - 
Chakravarthy- Osher limiter 
81x17 rnesh (48) 
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Figure 6: 

Pressure distribution 
Roee scheme, k = ^ 
Chakravarthy-Osher limiter 
81x17 mesh (48) 



Figure 7: 

Influence of the turbulence modeling 
Roe scheme 


a) Baldwin-Lomax model 


b) Johnson-King model 
321x65 mesh (192) 
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Figure 12: 

Pressure distribution (RAE2822) 
Roe scheme 
Harten-Osher limiter 
81x17 mesh (48) 




Figure 13: 

Pressure distribution (RAE2822) 
Roo scheme 

Chakravarthy-Oshcr limiter k = ^ 
81x17 mesh (48) 


Figure 14: 

Pressure distribution (RAE2822) 
Roe scheme 

Unlimited third order scheme 
81x17 mesh (48) 
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Figure 20: 

Pressure (a) and Skin friction (b) distributions 
Chakravartliy-Osher limiter k = ^ 
161x33 rnesli (96) 
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Figure 21: 

Pressure (a) and Skin friction (b) distributions 
Unlimited k = | 

161x33 mesh (96) 
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Figure 25: 

Pressure distribution 
Unlimited k = — 1 
161x33 mesh (96) 
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Figure 26: 

Convergence history 
Unlimited k = -1 
CFL=30, relax=0.5 
161x33 mesh (96) 
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Figure 32: 

Pressure distributions 

a) unlimited k = — 1 

b) “Smooth limiter” k = 5 

c) Chakravarthy-Osher limiter k = 5 

d) unlimited k = 3 

81x17 mesh (48) 
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Figure 33: 

Pressure and skin friction distributions 

a) unlimited k = — 1 

b) “Smooth limiter” k = 5 

c) Chakravartliy-Osher limiter k = ^ 

257x65 mcsh(192) 
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Figure 34: 

Pressure and skin friction distributions 
RAE2822 (case 10) 

“Smooth limiter” k = 3 
161x33 mesh (96) 



Figure 35: 
Convergence history 
Rae2S22 (case 10) 
“Smooth limiter” k = 
C’FL— 30, relax =0.5 
161x33 mesh (96) 
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Figure 36: 

Convergence history 

ONER A M6 Me,, = .84, a = 3.06, Re = 11 x 10 6 
a.) unlimited k = — 1 

b) “Smooth limiter” n — ^ 

c) Chakravarthy-Osher limiter k = k 

97x25x17 mesh 
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Figure 37: 

Pressure distributions 

ONER A M6 ilioo = .84, « = 3.06, Re = 11 x 10 6 
97x25x17 mesh 
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Figure 38: 

P rcssure dist ributi on s 

ONERA M6 Moo = .84, a = 3.06, Re = 1 1 x 10' 
Baldwin- Lomax model 
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Figure 39: 

Pressure and skin friction distributions 
Rae‘2822 (case 10) 
Baldwin-Lomax model 

a) 81x17 mesh (48) 

b) 161x33 mesh (96) 


c) 257x65 mesh (192) 






Figure 40: 

Pressure and skin friction distributions 
Rae2822 (case 10) 
Johnson-King model 

a) 161x33 mesh (96) 

b) 2. r )7x65 mesh (192) 

c) 321x65 mesh (256) 
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Figure 41: 

Pressure distribution 
Rae‘2822 (case 10) 
original Johnson-King model 
161x33 mesh (96) 
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Figure 42: 

Pressure distribution 
Rae2822 (case 10) 

J-K model with B-L formulation for v i% 
161x33 mesh (96) 



x/c 


Figure 43: 

Pressure distribution 
Rae2822 (case 10) 

J-K model with B-L formulation for v t 
257x65 mesh (192) 






Figure 44: 

Convergence history 
RAE2822 (case 10) 
Johnson-King model 
257x65 mesh (192) 
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Figure 45: 

Pressure and skin friction distributions 
RAE2822 (case 9) 
Johnson-King model 
257x65 mesh (192) 
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Figure 49: 

Pressure distributions 

ONERA M6 = .84. a = 6.06, Re = 11 x 10 6 
Baldwin- Lomax and Johnson-King models 
193x49x33 mesh 



Figure 50: 

Wall streanlines 

ONERA M6 = .84, o = 6.06, Re = 11 x 10 6 
Johnson-King model 
193x49x33 mesh 



Figure 51: 

Pressure contours 

ONERA M6 M 0 0 = .84, a = 6.06, Re = 11 x 10 6 
Johnson-King model 
193x49x33 mesh 
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Figure 53: 

Pressure contours 

ONERA M6 Moo = .84, a = 6.06, Re = 11 x 10 6 
Johnson-King model 
97x25x17 mesh 


Figure 54: 

Wall streamlines 

ONERA M6 Moo = .84, a = 6.06, Re = 11 x 10 6 
Jolinson-King model 
97x25x17 mesh 
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Figure 56: 

Pressure contours 

ONER A M6 Moo = .84, a = 3.06, Rc = 11 x 10 6 

a) Baldwin-Lomax model 

b) Jolmson-King model 

193x49x33 mesh 
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